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CHAPTER  1 
INTRODUCTION 


Loud  noises  such  as  gun  blast  noises  from  Army  training  facilities  have  brought 
concerns  from  nearby  communities.  In  response  to  this  environmental  issue,  the  Army  has 
built  and  implemented  noise  monitoring  systems  and  noise  barriers.  For  an  accurate  and 
effective  solution  to  noise  problems,  fundamental  acoustic  scattering  and  propagation  should 
be  carefully  studied. 

Studying  the  acoustic  scattering  has  been  difficult  and  incomplete  by  theoretical 
means.  An  efficient  and  accurate  numerical  method  is  necessary  to  investigate  acoustic 
scattering.  A  comprehensive  investigation  was  done  by  Alona  Boag  under  the  direction  of 
George  W.  Swenson,  Jr.  at  United  States  Army  Construction  Engineering  Research 
Laboratory  (USACERL)  Acoustics  Team.  Alona  Boag  concluded  that  method  of  moments 
(MOM)  was  the  best  approach  in  solving  acoustic  scattering  problems.  MOM  involves  in 
solving  an  integrodifferential  equation  in  a  matrix  equation  form.  She  has  written  an  MOM 
code  which  calculates  the  radiatioft  patterns  of  an  acoustic  source  scattered  by  an  object.  Her 
simulation  result  closely  matches  experimental  data  taken  from  a  scattering  of  a  rectangular 
baffle. 

Her  MOM  program  is  an  excellent  way  of  solving  scattering  problems  of  small  and 
simple  objects.  However  her  geometry  meshing  scheme,  explained  in  Chapter  2,  becomes 
tedious  when  creating  curved  or  complex  geometry.  To  effectively  create  3-D  objects,  a 
package  called  PATRAN  has  been  explored.  A  few  simple  PATRAN  commands  create 
geometry,  mesh  surfaces,  and  output  node  coordinates. 

Another  difficulty  is  in  analyzing  an  object  large  compared  to  its  wave  length,  X.  The 
number  of  unknowns  to  be  solved  grows  as  an  order  of  (d/X)^,  and  the  central  processing 
unit  (CPU)  time  grows  as  an  order  of  (d/X,)^,  where  d  represents  a  dimension  of  the  object. 
For  axisymmetric  geometry,  a  method  called,  body  of  revolution  (BOR)  technique. 
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eliminates  such  problems.  BOR  reduces  the  number  of  unknowns  to  an  order  of  (d/X), 
thereby  reducing  the  CPU  time.  A  BOR  code  written  in  FORTRAN  is  used  to  analyze 
scattering  from  curved,  axisymmetric  bodies.  The  BOR  code  has  been  extensively  tested  for 
its  validity,  including  a  comparison  to  an  exact  solution  for  scattering  from  oblate  spheroidal 
objects  involving  expansion  of  the  field  in  series  of  spheroidal  functions. 

Chapters  2  and  3  of  this  thesis  present  the  study  of  MOM  and  BOR  in  detail: 
mathematical  derivations,  code  implementation,  and  simulation  results.  In  addition,  an  in 
depth  study  of  oblate  spheroidal  functions  is  given  in  chapter  4.  Chapter  5  presents  the 
overall  results  and  conclusion. 
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CHAPTER  2 

METHOD  OF  MOMENTS 


There  are  several  methods  for  solving  an  acoustic  scattering  problem.  A  well  known 
method  is  geometrical  optics  (GO).  GO  yields  an  approximate  solution  for  scattering  from 
bodies  large  compared  to  the  acoustic  wavelength  [1].  Because  it  is  an  approximate  solution, 
only  applicable  to  the  analysis  of  high  frequency  waves,  it  is  inappropriate  for  USACERL’s 
scattering  analysis.  For  example,  a  typical  spectrum  of  an  Army  gun  blast  noise  ranges  from 
5  Hz  to  60  Hz  [2].  In  addition,  accurate  prediction  of  scattering  from  a  complicated  geometry 
is  almost  impossible.  The  above  reasons  eliminate  the  further  study  of  GO. 

Another  well  known  scattering  mode  is  Rayleigh  scattering.  Rayleigh  scattering  is 
typically  used  to  analyze  scattering  from  small  particles  such  as  air  bubbles  in  liquids  and 
water  molecules  in  the  air.  The  dimensions  of  practical  noise  barriers  are  comparable  to  or 
larger  than  the  acoustic  wavelength,  so  Rayleigh  scattering  is  not  an  appropriate  model. 

Alona  Boag  has  investigated  two  rigorous  numerical  methods  in  solving  acoustic 
scattering:  a  partial  differential  equation  method  (PDE)  and  an  integral  equation  method  (BE). 
She  presented  her  findings  at  a  USACERL  seminar,  and  Table  2.1  summarizes  her  findings. 
Table  2. 1  clearly  indicates  that  the  IE  method  is  the  best  choice  among  available  scattering 
solutions.  The  BE  method  leads  to  MOM,  and  the  remainder  of  this  chapter  is  devoted  to 
MOM. 

2.1  Mathematical  Derivation  of  The  IE 

When  an  acoustic  source  radiates  a  wave  in  free  space,  the  radiated  field  \|/  at  an 
observation  point  r  can  be  mathematically  described  as : 

\|r(r)  =  -  J  f(rs)  G(r,rs)  dv  (2.1) 

where  f  is  the  source  at  r*.  The  Green’s  function  G(r,rs)  is  an  impulse  response  at  r  due  to 
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f(rs).  In  a  homogeneous  medium. 


^gjkR) 

47tR 


(2.2) 


where  R  =  |r  -  FjI  .  The  wave  number  k  is  equal  to  2kIX.  The  integration  is  taken  over  the 
entire  volume  v  of  the  source.  Unless  otherwise  noted,  time  dependence  exp(-io)t)  is 
assumed. 

Table  2.1  Alona  Boag’s  comparison  between  PDE  and  IE 


Method 

Partial  Differential  Equation 

Integral  Equation 

Advantages 

•  applicable  to  bodies  comparable  to 
or  smaller  than  the  wavelength 

•  can  analyze  arbitrary  media  and 
geometry,  including 
inhomogeneous  bodies 

•  sparse  matrices  to  invert 

•  IE  is  based  on  the  Sommerfeld 
Radiation  Condition 

•  the  radiation  condition 
automatically  satisfied 

•  lower  number  of  unknowns 
compared  to  PDE  method 

•  high  accuracy 

Disadvantages 

•  number  of  unknown  N  grows  as 
(d/X)^ 

•  requires  artificial  boundary 
conditions  to  truncate  the  mesh 

•  medium  accuracy 

•  number  of  unknown  N  grows  as 
(d/X)^ 

•  applicable  only  to  piecewise 
homogeneous  problems 

•  CPU  time  grows  as  (d/X)^ 

In  the  presence  of  an  object  in  the  medium,  a  portion  of  the  radiated  wave  encounters 
the  object  and  re-radiates  into  different  directions.  The  object  acts  as  a  scatterer  of  the  original 
wave.  When  a  field  is  measured  at  r,  it  is  a  sum  of  the  fields  due  to  the  original  wave  and  the 
re-radiated  one.  The  field  due  to  the  original  wave  is  called  an  incident  field,  and  the  other 
field  due  to  the  re-radiated  one  is  called  a  scattered  field.  The  sum  of  the  incident  field  and  the 
scattered  is  the  total  field.  To  distinguish  each  field,  the  incident,  scattered,  and  total  fields  at 
r  are  denoted  as  P*,  and  P  respectively.  Thus  the  field  in  (2.1)  should  be  re-expressed 
as  the  incident  field  The  incident  field  is  recalculated  using  the  Green’s  function 
assuming  that  the  presence  of  the  scatterer  does  not  influence  the  propagation  of  the  original 
wave. 
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The  homogeneous  Green’s  function  satisfies  the  Helmholtz  equation,  i.e.,  the 
frequency  domain  wave  equation: 

jvS^)\K(r)  =  f(rs)  (2-3) 

The  speed  of  the  sound  is  denoted  as  c  and  is  assumed  to  be  340  m/sec  in  the  air. 
Substituting  (2.1)  in  (2.3)  yields, 

(v"+^|G(r,rs)  =  -8(r-r,)  (2.4) 

Multiplying  (2.3)  by  G<r,r,)  and  (2.4)  by  tKr,r,)  and  rearranging  terms  in  each  equation, 
the  following  relations  are  obtained: 

G(r,rs)  V  V(r)  =  G(r,rs)  jf(rs)  -  ^  \l/(r)J 

\|/(r)  V^G(r,rs)  =  -  V(r)  [sfr-Fj)  +  ^  G(r,rs)j 
The  Green's  second  identity  states  that  for  arbitrary  scalars  U  and  V, 

£(uvV-vvMdn=jju|^.v§)ds  (2.7) 


(2.5) 

(2.6) 


where  Q  is  a  volume  and  S  is  the  surface  enclosing  Q.  Replacing  U  with  the  pressure  field, 
V  with  the  Green’s  function  and  using  relations  (2.5)  and  (2.6), 


5(r-rs)  G(r,rs)  +  ^  G(r,rs)]  +  G(r,rs)  (f(rs)  -  ^  W) 
c2  /  V  C^ 


dQ  . 


3G(r,r') 

9n 


9tKr’) 
G(r,r’)  ^ 


9n 


dS 


(2.8) 


where  \ir(r')  is  the  field  on  the  surface  and  the  partial  derivatives  were  taken  with  respect  to  a 
normal  vector  on  the  surface.  The  normal  vector  n  points  outward  from  the  volume  Q.  Note 
that  the  minus  sign  at  the  right  hand  side  of  (2.7)  reflects  that  -n  is  taken  instead  of  n  (see 
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Figure  2.1).  The  surface  includes  the  scatterer  surface  and  the  fictitious  surface  at  infinity 
(see  Figure  2.1). 


Figure  2.1  A  scatterer  and  a  source  in  an  unbounded  space 

After  substituting  (2.5)  and  (2.6)  on  the  left  hand  side  of  (2.8)  and  simplifying,  (2.8) 
becomes, 

-V(r)  -  f  G(r,rs)  f(rs)  dfl  =  -  f  f -\|/(r')  -  G(r,r')  dS  (2.9) 

Js  V  dn  j 

By  equation  (2.1),  the  integral  on  the  left  hand  side  of  (2.9)  is  essentially  equal  to  the 
incident  field 

All  fields  in  any  unbounded  medium  must  satisfy  the  Sommerfeld  radiation  condition; 

a  \ 

- ik  \|/ =  0  as  R ->  oo  (2.10) 

9R 

This  condition  ensures  that  the  field  \|r  becomes  negligible  at  a  large  distance  from  the 
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acoustic  source.  Applying  the  radiation  condition  to  (2.10)  eliminates  the  field  contribution 
from  the  fictitious  surface.  Thus  the  volume  integral  is  reduced  to  a  surface  integral: 


P(r)  =  pinc(r)  +  p(r') 


aG(r,r') 

3n 


G(r,r') 


3P(r') 


0n 


dSs 


(2.11) 


where  \|/(r)  and  \|r(r')  are  re-expressed  as  the  pressure  field  P(r)  at  r  and  P(r')  at  r’ 
respectively. 

P(r),  P(r'),  and  the  partial  derivative  of  P(r')  are  the  three  unknowns  to  be  solved. 
The  surface  is  assumed  to  be  perfectly  rigid.  Rigid  surfaces  do  not  move  in  the  presence  of 
acoustic  pressure,  thus  the  normal  velocity  field  on  the  surfaces  is  zero.  The  pressure  field  on 
the  surface  is  doubled,  and  its  partial  derivative  normal  to  the  surface  is  zero.  Finally  the  IE 

of  interest  is 

P(r)  =  F"‘^(r)  + 


1 


P(r')^^^dSs 
0n 


(2.12) 


(2. 12)  will  be  used  to  build  a  matrix  equation  that  can  be  solved  in  a  computer  program. 


2.2  Evaluation  Of  The  Singularity 

The  equation  (2.12)  can  be  evaluated  at  any  observation  point  r,  and  when  r 
approaches  the  surface  of  the  scatterer,  the  integral  in  (2.12)  becomes  singular  at  r— r . 
Morita,  et.  al.  cleverly  evaluate  the  integral  at  the  singularity  [3]. 

Take  an  infinitesimally  small  surface  around  r'  and  call  it  5S.  For  a  convenient 
evaluation  of  the  singularity,  shift  the  coordinate  system  such  that  the  observation  point  r  and 
the  integration  point  r'  are  (0,0,z)  and  (p’,(t)’,z'),  respectively  (see  Figure  2.2)  and  set  z’=0. 
The  integral  is  then  approximated  as. 
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At  r=r',  P(r')  is  approximately  equal  to  P(r)  and  may  be  assumed  invariant  within  5S. 
Figure  2.2  indicates  that  -z  is  the  normal  vector,  and  taking  a  derivative  with  respect  to  z 
yields, 


1 

+  (z-z')^ 


2  TCp' 


dp'  = 


(2.14) 


r 


Figure  2.2  Limiting  procedure 

As  a  and  z  approach  zero,  the  integral  at  the  singularity  is  reduced  to  -  ^  P(r).  Therefore 

when  the  observation  point  is  taken  on  the  surface  of  the  scatterer,  (2.12)  should  be  rewritten 
as: 

pinc(r)=Ip(r)+j'  p(r-)^G^dSs  (2.15) 

where  f  indicates  the  principal  value  of  the  integral.  The  above  equation  is  solved  to  find  the 

pressure  fields  on  the  surface  by  assembling  (2.15)  into  a  matrix  equation.  Once  the  pressure 
fields  on  the  surface  are  known,  then  the  real  IE  (2.12)  can  be  solved  for  any  observation 
location. 


2.3  Matrix  Equation  Assembly 

Assembling  a  matrix  equation  from  IE  is  the  key  step  in  MOM.  MOM  is  a  technique 


which  solves  the  EE  by  dividing  the  scatterer  surface  into  small  patches  called  elements.  The 
elements  can  be  any  shape  as  long  as  the  shape  can  be  used  tightly  to  mesh  the  surface, 
leaving  no  “holes”  or  space  on  it.  The  simulation  program  written  by  Alona  Boag  uses 
triangular  elements  to  mesh  the  surface.  For  accurate  simulation  results,  the  lengths  of  three 
sides  of  each  triangular  element  must  be  approximately  equal,  i.e.,  isometric.  Creating  an 
element  with  one  extremely  long  leg  compared  to  the  others  is  not  recommended.  Also  the 
length  of  each  triangle  leg,  denoted  by  grid  length,  should  be  less  than  one  tenth  of  a  wave 

length. 

Assuming  the  pressure  fields  do  not  very  drastically  within  each  element,  P(r )  in 
(2. 15)  can  be  regarded  constant  on  each  element.  Usually  the  value  of  P(r )  on  the  center  of 
each  element  is  picked;  this  is  how  the  method  of  moment  is  named.  The  process  of 
assembling  a  matrix  equation  is  described  for  a  typical  rectangular  barrier  (see  Figure  2.3). 


z 


Figure  2.3  The  meshed  rectangular  barrier  and  the  coordinate  system 


Let  N  be  the  total  number  of  the  triangular  elements  on  the  surface,  and  define  pulse  basis 


functions 


fj  =1,  on  AS’j 

fj  =  0,  on  all  other  AS’i 


(2.16) 


where  AS'i  is  the  j**’  triangular  element  on  the  surface.  Let  r'j  be  the  center  coordinate  of 
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AS'i  ■  Also  let  the  unknown  pressure  fields  on  the  surface  be  expressed  by, 


N 


P(r')  =  I  Pj  fj 

i=i 


(2.17) 


where  p  is  a  1  x  N  unknown  array  for  the  surface  pressure  fields.  Substituting  (2.17)  into 
(2.15),  and  matching  the  field  at  the  midpoint  of  each  A’s,  a  matrix  equation  is  obtained  as 
below: 


N 


Pl"^  =  I  Zij  Pj 
i=l 

Expressing  (2.18)  in  a  matrix  form. 


,  i=l,2,3,..,N 


(2.18) 


pine 

(2.19) 


where 


0G(ri,r') 

0n 


dS 


1 

2 


i=j 


(2.20) 


The  derivative  of  the  Green’s  function  with  respect  to  a  normal  vector  is 

.  =  (ik _ 1  \exp(ik|ri- r'l) 

r  Iri-r’l)  47r|r|-rl 


aG 

an 


(2.21) 


The  N  X  N  matrix  Z  is  called  an  impedance  matrix,  and  it  needs  be  inverted  to  find  the 
unknown  surface  pressures  in  p. 


N 


PJ  =  I  Zij  ri"' 


i=l 


(2.22) 


To  find  a  total  pressure  at  any  location  r,  sum  the  incident  field  and  the  scattered  field  found 
by  the  above  method; 
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(2.23) 


N  f 

P(r)  =  P'"‘'(r)  +  I  pj 

j=l  JASj 


30^  ds 

3n 


2.4  Impedance  Matrix  Formulation 

Finding  values  of  the  Z  matrix  is  the  central  part  for  solving  the  matrix  equation. 
Each  matrix  element  Zij  represents  the  contribution  of  the  surface  integration  of  element  at 
the  center  of  i*  element.  The  surface  integration  is  numerically  done,  and  its  accuracy  is  very 
critical  to  the  overall  accuracy  of  MOM.  The  seven-point  Gaussian  quadrature  is  used  for  the 
numerical  integration;  an  integration  of  any  ‘smooth’  function  within  a  triangle  can  be 
accurately  obtained  by  sampling  seven  points  within  the  triangle. 

f  f(r)  ds  =  A  X  Wi  f(ri)  +  O(h^)  (2.24) 

Jtriangle 

where  A  is  the  area  of  the  equilateral  triangle  and  wj  is  the  weighting  value.  The  error  of  the 
quadrature  is  on  the  sixth  order  of  h,  the  radius  of  a  circle  circumscribing  the  triangle,  so  h 
must  be  less  than  1  for  accuracy.  A  point  P  inside  a  triangle  can  be  uniquely  expressed  in 
terms  of  the  areas  of  subtriangles  defined  by  P  and  the  vertices  of  the  main  triangle  (see 
Figure  2.4).  The  subtriangles  LI,  L2,  and  L3  are  called  area  coordinates  and  they  add  up  to 
1.  The  area  coordinates  and  weights  for  the  seven  points  are  listed  in  Table  2.2  [4].  The 
coordinates  for  the  point  inside  the  triangle  can  be  expressed  as 

Ti  =  P 1  LI  +  P2  L2  +  P3  L3  (2.25) 


where  PI,  P2,  P3  are  the  three  vertices  coordinates  of  the  triangle. 


i 

LI 

L2 

L3 

Wj 

1 

1 

1 

1 

270 

3 

. .  3 

.  3 

HHBwiiiBlii 

2 

(9  +  2^) 

(4-VTy) 

(T-VTJ) 

(155 -flS) 

21 

21 

21 

1 200 

3 

(4-VT5) 

(9  +  2VT5) 

(7-VT5) 

(155 -■/T3') 

21 

21 

91 

i2nn 

4 

(4-VT5’) 

(4-/TJ) 

(13  +  2103') 

(155 -lOT) 

21 

21 

21 

1200 

5 

(9-2VT5') 

(4  +  YT5) 

(7  +  m') 

(155+fI3') 

21 

21 

21 

1200 

6 

(4  +  ^) 

(9  -  2fB) 

(7+VT3') 

(155+VT3') 

21 

21 

91 

1900 

7 

(4  +  iTS) 

(4  +  m) 

(13-2V1J) 

(155+101) 

91 

91 

91 

1200 

2.5  Simulation  Results 

To  verify  the  accuracy  of  the  MOM  program,  a  simulation  of  scattering  from  a 
wooden  baffle  was  run  and  compared  to  the  physical  measurements  performed  outdoors. 
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The  baffle  was  61  cm  wide  and  30.5  cm  high,  and  a  microphone  was  placed  at  the  horizontal 
center,  about  7.6  cm  above  the  ground,  and  2.5  cm  away  from  the  face  of  the  baffle.  The 
experimental  data  are  taken  from  Benson  et  al.  [5].  In  the  numerical  analysis  a  point  source  is 
placed  at  the  microphone  position  using  the  reciprocity  theorem.  The  simulation  results 
closely  match  the  measurements  (see  Figures  2.5  and  2.6).  Also  diffraction  patterns  for  the 
same  baffle  but  mounted  at  a  tilt  angle  have  been  simulated  (see  Figures  2.7  and  2.8). 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  2.5  The  MOM  baffle  simulation  and  measurements  1 
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Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  2.6  The  MOM  baffle  simulation  and  measurements  2 


Side  view 


+  Z-axis 


Figure  2.7  The  tilted  baffle 
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Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  2.8  MOM  simulations  of  the  tilted  baffle  2 


When  one  fires  a  gun  inside  a  shed,  the  acoustic  excitation  may  be  represented  as  a 
point  source.  A  simple  shed  has  been  modeled  as  in  Figure  2.9,  and  the  diffraction  patterns 
have  been  studied.  Figure  2.10  shows  the  diffraction  pattern  of  a  point  source  inside  a  shed. 
Also  a  simple  hill  has  been  modeled  to  see  the  diffraction  patterns  of  gun  blast  noises  nearby 
a  hill  (see  Figures  2.1 1  and  2.12).  The  dimensions  of  the  hill  is  small  compared  to  a  realistic 
hill,  and  analyzing  a  larger  hill  requires  more  CPU  time  and  memory. 
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Top  view 


1.5  m 


Side  view 


Figure  2.9  The  shed 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Top  view 


Front  view 
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Total  pressure  field  amplitude  in  the  azimuthal  plane 


- Point  Source  Position:  (0,  5.5,  0.5)  m  I 

- Point  Source  Position:  (0,  8.0,  0.5)  m  I 

Figure  2.12  The  MOM  simulations  of  the  hill 

2.6  Problems  In  MOM 

Although  Alona  Boag’s  program  accurately  predicts  the  scattering  pattern,  there  are 
several  difficulties  in  using  her  codes.  First  of  all,  the  mesh  generator  she  used  is  not  user- 
friendly.  For  example,  to  mesh  the  rectangular  baffle,  one  has  to  enter  four  vertices 
coordinates  for  each  rectangular  face  of  the  baffle  counter  clockwise,  looking  from  the 
outside  of  the  baffle.  This  ensures  the  correct  calculation  of  outward  normal  vectors.  It  is 
simple  for  the  baffle  with  six  faces,  but  when  a  geometry  is  complex  and  has  many  facets,  it 
is  tedious  to  enter  the  coordinates  manually.  In  addition,  an  accurate  meshing  of  any  curved 
surface  is  impossible  since  the  mesh  generator  assumes  a  flat  surface.  To  create  geometry 
and  mesh  its  surfaces  with  ease,  a  finite  element  analysis  package  called  PATRAN  is 
explored.  A  few  simple  PATRAN  commands  create  geometry,  mesh  surfaces,  and  output  the 


coordinates.  Creating  and  meshing  curvatures  are  simple  with  PATRAN.  However,  there  is 
a  drawback  in  PATRAN.  To  define  outward  normal  vectors  on  surfaces  of  an  object,  the 
coordinates  should  be  listed  in  counter  clockwise  sequence.  PATRAN  meshes  and  lists 
element  coordinates  arbitrarily,  either  clockwise  or  counter  clockwise.  Some  normal  vectors 
point  outward,  and  others  inward.  To  remove  this  arbitrariness,  a  simple  routine  is  written. 
Assuming  each  element  is  approximately  equal  in  size,  the  center  coordinates  of  each  elemait 
are  summed  up  and  then  divided  by  the  number  of  elements,  N.  This  yields  an  approximate 
center  coordinate  of  the  scatterer,  and  a  dot  product  of  a  vector  from  this  center  of  the  object 
to  a  center  of  any  element  and  the  corresponding  normal  vector  should  be  a  positive  quantity. 
If  negative,  then  it  indicates  the  normal  vector  is  pointing  toward  the  center  of  the  object,  and 

the  direction  of  the  vector  should  be  reversed. 

Other  difficulties  are  the  memory  and  CPU  time  needed  in  computing.  Since  the 
number  of  elements  N  is  proportional  to  the  surface  area,  analysis  of  an  object  large 
compared  to  the  wavelength  is  difficult,  if  not  impossible.  As  N  increases  the  required 
memory  for  storing  the  impedance  matrix  increases  as  N^.  In  FORTRAN,  reserving  a  large, 
contiguous  section  of  working  memory  becomes  difficult  when  solving  for  a  large  number  of 
unknowns.  For  example,  in  the  Michigan  baffle  simulation  at  163  Hz,  7,308  triangle 
elements  are  meshed,  and  the  number  of  the  matrix  elements  is  53,406,864.  Each  matrix 
element  occupies  16  bytes  of  memory  for  a  complex  double  precision  data  type.  Thus  more 
than  854  MBytes  are  needed  in  creating  the  matrix.  The  CPU  time  for  a  matrix  inversion  is 
proportional  to  N^  Most  of  the  acoustic  scattering  simulations  have  been  solved  in  a  Convex 
C420  machine  at  the  University  of  Blinois.  The  Convex  C420  has  four  processors  with  512 
MBytes  of  memory.  It  is  used  for  high  speed,  heavily  vectorized  computations.  However 
C420  has  proved  to  be  inadequate  for  the  baffle  simulation  at  163  Hz.  An  access  to  a  mofe 
powerful  machine  is  obtained  from  National  Center  for  Supercomputing  Applications 
(NCSA):  Power  Challenge  account.  The  Power  Challenge  Machine  is  built  by  Silicon 
Graphics,  Inc.,  and  it  has  16  shared  memory  multiprocessors.  Its  total  memory  is  four  giga 
bytes,  and  it  can  perform  300  mega  floating  point  operations  per  second  (flops)  per 
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processor,  total  of  4.8  giga  flops. 

In  spite  of  the  increased  capacity,  a  more  powerful  machine  is  needed  in  analyzing 
larger  objects  at  high  frequencies.  The  limited  resources  are  unavoidable.  One  should  look 
for  a  clever  method  of  solving  the  problem  instead  of  looking  for  a  faster  machine.  A 
technique  that  overcomes  all  of  the  above  problems,  meshing,  memory,  and  CPU  time,  is 
sought  and  implemented  in  Chapter  3. 
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CHAPTER  3 

BODY  OF  REVOLUTION 


As  discussed  in  Chapter  2,  MOM  has  many  limitations  in  analyzing  large,  curved 
scatterers.  A  common  example  of  such  a  scatterer  is  a  parabolic  reflector.  A  parabolic 
reflector  is  axisymmetric,  so  a  method  called  body  of  revolution  (BOR)  can  be  used  to  solve 
such  a  case.  The  study  of  the  BOR  technique  is  based  on  a  paper  by  Seybert  et  al.  [6].  By 
taking  an  advantage  of  the  axisymmetric  properties,  the  surface  integral  is  reduced  to  a  line 
integral  along  the  generator  of  the  scatterer  body  and  an  integral  over  the  angle  of  the 
revolution.  The  integration  over  the  angle  is  performed  partly  analytically  in  terms  of  elliptic 
integrals  and  partly  numerically  using  the  Gaussian  quadrature  formula.  A  program  in 
FORTRAN  is  written,  and  the  same  test  cases  in  the  Seybert’ s  paper  are  run  and  compared  to 
the  paper.  My  results  agree  well  with  both  the  paper  and  the  theoretical  solution.  This  chapter 
describes  BOR  in  detail. 


3.1  Derivation  Of  BOR 

The  derivation  of  BOR  is  based  on  the  Helmholtz  integral  formula. 


C(r)  \i<r)  =  \|/(r') 


(3.1) 


where 


C(r)  =  47t  + 


rii 

1  \ 

Is^ 

lR(r,r')) 

ds(r’) 


(3.2) 


Throughout  the  derivation  of  BOR,  the  time  dependence  of  exp(i(Ot)  is  used,  and  the  Green’s 
function  is  taken  as  G(r,r')  —  exp(-ikR)/R  to  follow  the  notations  of  the  Seybert  paper.  The 
integral  in  (3.2)  approaches  -27t  on  the  smooth  surfaee  of  the  body. 

Figure  3.1  illustrates  a  simple  axisymmetric  scatterer  and  the  nomenclature  needed  in 
the  derivation  of  BOR.  When  the  scatterer,  boundary  conditions,  and  sources  are 
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axisymmetric,  the  pressure  fields  and  the  normal  velocity  fields  on  the  surface  of  the  scatterer 
are  also  axisymmetric.  Taking  an  advantage  of  this  symmetric  distribution  of  the  fields,  the 
surface  integral  may  be  reduced  to  a  line  integral. 
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The  above  equation  can  be  decomposed  into  two  parts 
Ki(r 


_  f  '^_^exp(-ikR(r.r’)  4|  Q 
M  R(r,r')  ' 


(3.5) 


and 


K2(r,r') 


)»27t 
0  ' 


1 


9nlR(r,r')| 


k\  <i9(f) 


(3.6) 


As  R  approaches  zero,  the  integrand  of  Ki(r,r')  can  be  expressed  in  terms  of  a  Taylor’s 
series  expansion, 

exp(-ikR(r,r')  -  1  _  1  -  ikR(r,r')  +  0.5  (-ikR(r,r'))2  -  1  ^  R(r,r')  ^3 

R(r,r')  R(r,r')  2 


Let  the  normal  component  of  R(r,r')  be  d,  and  the  tangential  component  be  x.  Then  the 
derivative  of  (3.7)  becomes, 

8[.  k2R(r,rM_  3R(r,r')  _  d  (3.8) 

■  2  /  ■  2  9n  2  Vx2  +  d2 

As  X  approaches  zero,  (3.8)  will  be  so  Ki(r,r')  is  non-singular  and  can  be  evaluated 
numerically  at  R(r,r')  =  0.  K2(r,r')  can  be  expressed  in  terms  of  the  elliptic  integrals. 


K2(r,r-)  =  ^ 


rln  ^ 

Jo  R(r,r') 


d0(r') 


a_ 

dn 


^2k 

.0 


dtp 


i?  sin^tp) 


1/2 


8n 


4F(7i/2,k) 

R 


(3.9) 


where  F(7t/2,  k)  is  a  complete  elliptic  integral  of  the  first  kind  with  modulus  k .  Also, 

^  _  (Pq  +  pp)2  4.  (zq  +  zpf  and  k  =  The  outward  normal  vector  at  point  Q  is 

R 


defined  as  n  =  npQp  +  nzQ^,  and  the  angle  of  revolution  is  0  =  0q  -  0p.  Evaluating  (3.9) 
yields. 


K2(r,r')=:^ 

r2 


U9F  8k  p  aRL  . 
^^8k^pQ’  ^Pq) 


=-8F  8k  p  8R\ 

i  8kazQ'  8zq) 


(3.10) 


The  derivative  of  F  is  expressed  in  terms  of  the  complete  elliptic  integral  of  the  second  kind. 
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(3.11) 


ap 


_  E  -  k^F 


ak  kk^ 

where  E  =  I  (l  -  ^sin^cp)  d(p  and  k’  =  (l  - Substituting  above  expressions  and 

Jo 

evaluating  derivatives  in  (3.10), 


K2(r,r')  =  ^ 


^1  RE  1  RF  E(PQ  +  Pq)]j^  E(ZQ-Zp) 
2pqF  2  Pq  rF  j  PQ  Pr 


Thus  the  surface  integral  in  (3.1)  is  expressed  as  a  line  integral  as  below: 

C(r)  \j<r)  =  \|/(r')  K(r,r’)  p(r')  dL(r')  +  47t  \|i*"‘^(r) 

where  C(r)  =  K2(r,r')  p(r’)  dL(r')- 


(3.12) 


(3.13) 


3.2  Numerical  Implementation 

Equation  (3.13)  can  be  numerically  evaluated  along  a  contour  L  of  an  axisymmetric 
body.  The  line  L  is  divided  into  segments  for  integration  (see  Figure  3.2). 


z  Element  1 


Figure  3.2  Discretization  on  the  contour  L 
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The  line  L  is  divided  into  four  elements  where  each  element  is  made  up  of  three  nodes.  The 
contour  along  the  line  L  is  assumed  to  be  varying  quadratically,  and  the  variation  is  called  a 
shape  function,  Na(^).  A  quadratic  shape  function  is  defined  as  below,  and  it  provides  a 
more  accurate  description  of  the  contour  than  does  a  linear  shape  function.  ^  is  called  a  local 
coordinate  and  ranges  from  -1  to  1,  and  a  is  called  a  local  node  number  and  ranges  from  1  to 

3. 


,  a=  1 

N„{4)  =  1  - 

,  a  =  2 

(3.14) 

,  a  =  3 

The  actual  coordinates  can  be  represented  in  terms  of  the  shape  function  and  local  coordinate, 

p(^)  =  X  Na(^)  pa 

(3.15) 

z(^)=  i  Na(^)Za 

a=l 


The  above  equations  map  each  element  on  the  body  onto  a  straight  line  in  the  local  coordinate 


Element  i  Local  coordinate 

Figure  3.3  The  element  mapping  onto  the  local  coordinate  system 

Similarly  the  pressure  fields  \|r  along  L  are  approximated  using  the  same  shape  function. 

Vm(^)  =  X  Na(4)  Vma 

«=i  (3.16) 

where  \|rma  is  the  value  at  node  a  on  the  element  m.  Because  the  fields  distribution  and  the 
contour  are  described  using  the  same  shape  function,  the  formulation  is  said  to  be 
isoparametric. 


25 


Applying  the  above  scheme  to  the  line  integral  Equation  (3.13)  yields, 


N  3 


■i: 


C(r)iKr)=I  Iij/mal  K(r,r'(4))  Na(^)  p(^)  Jm(4)  V"=(r)  (3.17) 

m=l  a=l 


where 


3 

n=l  J-i 


C(r)  =  4jt  +  X  K2(r,r'(^))  p(^)  Jm(^)  d^ 


(3.18) 


N  is  the  total  number  of  elements  used  to  discretize  the  scatterer  body,  and  Jm(^)  is  the 
Jacobian  of  transformation  for  element  m. 

Jm  a)  = 


1/2 


^(p3  -pl)^  +  (pl  +p3  -2p2)2  +  i(z3-zl)^  +  (zl  +z3-2  z2)2 


(3.19) 


Therefore  for  any  point  P  on  the  surface  of  the  scatterer,  whose  global  node  number  is  j,  the 
integral  equation  can  be  expressed  as  follow: 


where 


and 


-i  i  Vmaa“j  +  Vjfl+  i  Cj 

(3.20) 

m=l  a=l 

\  m=l  f 

K(Pj,4)N„(^)  p(^)Ja^)d^ 

(3.21) 

c  ■  -  — 

4;: 

j  K2(Pj,4)  p(^)  J^(^)  d^ 

(3.22) 

Equations  (3.20)  to  (3.22)  apply  to  each  node,  and  those  M  equations  make  up  a 
matrix  equation  for  BOR. 


A 

= 

ync 

. 

(3.23) 


For  integration  on  the  ^  axis,  the  Gaussian  quadrature  was  used.  Initially  a  4-point  Gaussian 
quadrature  formula  was  used  to  evaluate  the  integral.  Comparing  with  the  theoretical  solution 
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revealed  that  the  4-point  Gaussian  quadrature  introduced  20  -  30%  error.  6-  or  8-point 
Gaussian  quadrature  must  be  used  for  an  accurate  numerical  integration.  In  addition  the 
number  of  total  nodes  were  varied,  and  the  number  of  nodes  did  not  contribute  much  to  the 
accurate  analysis.  For  example,  scattered  patterns  from  9,  17,  and  33  node-bodies  were 
similar  to  one  another.  Thus  beyond  a  reasonable  number  of  nodes,  the  number  of  total 
nodes  will  lengthen  the  computational  time  without  any  improvement  on  the  accuracy. 


3.3  Scattering  From  Spheres 

Several  cases  were  run  to  verify  the  validity  of  the  code.  To  compare  with  the  results 
from  the  reference  paper  [6],  identical  test  parameters  were  chosen.  Seybert  et  al  presented 
the  scattered  pressure  field  patterns  of  spheres  due  to  an  incident  plane  wave  [6].  The  test 
case  results  closely  agree  with  the  results  from  the  paper.  In  addition,  a  theoretical  scattering 
pattern  for  Figure  3.4  was  computed  in  Mathematica  for  a  further  comparison.  Due  to  its 
simple  geometry,  scattering  by  a  sphere  can  be  solved  analytically.  For  a  rigid  sphere  with  a 
plane  wave  incident  in  -z  axis,  the  total  field  is  [7], 


yot  =  +  \|/*  =  X  (-i)"  (2n  +  1)  |jn(kr)  -  a'n  h^n  ^(kr)]  Pn(cos0)  (3.24) 

n=0 


where  jn  and  h|/^  are  spherical  Bessel  and  Hankel  functions,  respectively.  Also  a'n 


j'n(ka) 

hj,^>'(ka) 


where  a  is  the  radius  of  the  sphere.  Pn  is  the  associated  Legendre  function.Since  the  incident 
field  is  exp(-ikr  cosO),  the  scattered  field  is  -  exp(-ikr  cos0).  In  the  actual  calculation  of 
the  analytic  solution,  the  summation  in  (3.24)  were  added  up  to  n  =  30  (see  Figure  3.5).  The 
test  cases  involving  only  one  scattering  body  were  compared  to  the  analytic  solutions,  and 
both  my  results  and  [6]  closely  followed  the  analytic  solutions  (see  Figure  3.6).  For  test 
cases  involving  two  discrete  spheres,  approximate  solutions,  which  neglected  the  interactions 
between  the  two  bodies,  were  compared  to  the  numerical  solutions  (see  Figures  3.7,  3.8, 
and  3.9). 
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Incident  planewave 


Figure  3.4  Scattering  from  one  sphere 


Scattered  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.5  The  BOR  simulation  and  analytic  scattering  patterns  of  the  sphere 


Scattered  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.6  The  BOR  simulation  of  the  sphere 


t  Incident 
planewave 
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Scattered  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.8  The  BOR  simulations  of  the  two  spheres  1 


Scattered  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.9  The  BOR  simulations  of  the  two  spheres  2 


3.4  Scattering  From  Parabolic  Reflectors 

Comparing  the  simulation  results  in  Section  3.3  is  not  adequate  for  a  complete 
validation  of  the  BOR  code.  Scattering  from  the  spheres  is  axisymmetric  regardless  of  the 
propagation  direction  of  the  incident  field.  Also  only  the  diffractions  due  to  a  plane  wave 
incidence  have  been  compared  and  validated.  For  a  complete  testing  of  the  BOR  code, 
scattering  due  to  a  point  source  incident  on  an  object  such  as  a  parabolic  reflector  should  be 

tested. 

Simulations  on  the  diffractions  of  the  parabolic  reflectors  can  be  useful  since  the 
parabolic  reflectors  are  used  as  acoustic  amplifiers  [8].  Recording  bird  sounds  in  the 
wilderness  is  not  a  trivial  task  for  ornithologists  since  the  bird  sounds  may  not  be  loud 
enough  to  be  distinguished  in  the  presence  of  background  noises.  Many  ornithologists  use 
parabolic  microphones  to  record  bird  sounds.  A  parabolic  microphone  is  a  parabolic  reflector 
dish  with  a  microphone  attached  to  the  dish.  Typically  the  microphone  is  located  at  the  focal 
point  of  the  parabolic  reflector  because  the  fields  will  arrive  in  equi-phase  at  the  surface  of  the 
reflector.  Using  the  reciprocity  theorem,  a  diffraction  pattern  due  to  a  point  source  at  a  focal 
point  is  identical  to  the  reception  pattern  of  the  microphone  located  at  the  focal  point. 

Although  the  use  of  the  parabolic  nucrophones  are  common  in  bird  recordings,  exact 
reception  patterns  are  not  supplied  by  their  manufacturers.  The  receptivity  of  parabolic 
reflectors  is  thought  to  be  highly  directional.  Manufacturing  parabolic  microphones  is  not  as 
simple  as  building  flat  baffle-type  directional  microphones  because  of  the  curvature.  Thus  the 
price  of  the  parabolic  microphones  are  quite  expensive.  The  diffraction  patterns  from  various 
parabolic  reflectors  can  be  accurately  calculated  by  the  BOR  code.  Note  that  as  the  focal  point 
is  place  farther  away  from  the  reflector,  the  curvature  of  the  reflector  becomes  more  flat  (see 
Figure  3.10).  Figures  3.11  through  3.15  show  far  field  diffraction  patterns  for  parabolic 
reflectors  whose  diameters  are  twice  the  acoustic  wavelength  of  the  interest.  The  reflector 
analyzed  in  Figure  3.15  is  virtually  a  disk  since  the  focal  point  is  100  wavelengths  away 
from  the  center  of  the  reflector.  Except  for  Figure  3.11,  the  highest  directivity  is  achieved 
when  the  point  source  is  placed  at  0.5  wavelength  away  from  the  center.  This  is  plausible 
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because  the  scattered  fields  on  the  surface  of  the  reflectors  have  the  same  phase  as  the 
incident  fields  on  the  rigid  reflectors.  Therefore  when  the  source  is  positioned  at  0.5 
wavelength  away  from  the  center,  the  incident  and  reflected  fields  are  also  in-phase,  adding 
up  the  field  amplitudes.  The  patterns  show  that  the  parabolic  reflectors  do  not  neccessarily 
achieve  high  directivity.  Large  side  lobes  can  be  found  depending  on  the  positions  of  the 
source.  Figure  3.16  compares  the  diffraction  patterns  of  a  parabolic  reflector,  disk,  and  a 
rectangular  baffle  with  similar  dimensions.  The  parabolic  reflector  achieves  slightly  higher 
directivity,  and  it  is  not  enough  to  justify  for  purchasing  expensive  parabolic  reflector 
microphones. 
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Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.11  The  BOR  simulations  of  the  parabolic  reflector  1 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.12  The  BOR  simulations  of  the  parabolic  reflector  2 
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Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.13  The  BOR  simulations  of  the  parabolic  reflector  3 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  3.14  The  BOR  simulations  of  the  parabolic  reflector  4 
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Figure* 


Suiirv  u  al  ihc'  suilaca 
•  Source  at  0.25  \Mwelength 
■  Source  al  0.5  wavelength 


Diameter:  2  wavelengths 

Focus:  100  wavelengths 
’’flat  disk” 


;  of  the  parabolic  reflector  5 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Parabolic  Reflectc*’ 

Diameter:  1  wavelength 
Focal  length:  0.22  wavelength 
Point  source  at  focus 


Disk 

Diameter:  1  wavelength 

Point  source  centM“ed  on  the  surface 


Square  plate 
Side:  1  wavelength 
Point  source  centered  on  the  surface 


Figure  3.16  The  BOR  simulations  of  the  parabolic  reflector,  disk,  and  baffle 


3.5  Scattering  From  Thin  Disks 

Although  the  parabolic  reflector  simulations  due  to  a  point  source  could  not  be 
compared  to  known  results,  the  simulations  can  still  be  checked  for  the  validity.  When  the 
focal  length  of  the  parabolic  reflector  becomes  large  compared  to  the  wavelength,  the  shape 
of  the  reflector  resembles  a  thin  disk.  Diffraction  of  a  thin  disk  has  been  rigorously  analyzed 
previously  by  many  scholars  [7,  9,  and  10].  Leitner  presented  the  diffraction  of  sound  by  a 
circular  disk  due  to  a  plane  wave  incidence  [7],  and  the  exact  theoretical  values  based  on  the 
wave  functions  of  the  oblate  spheroid  were  plotted  in  his  paper.  Also  Wiener’s 
measurements  supported  Leitner’ s  work  [10]. 

In  order  to  compare  with  the  Leitner’ s  results,  the  same  test  cases  were  run  using 
BOR  code  (see  Figure  3.17).  Figures  3.18,  3.19,  3.20,  3.21,  and  3.22  show  the  BOR 
simulation  results.  The  BOR  simulation  results  agree  with  the  Leitner’s  exact  theoretical 
calculations,  and  they  match  closer  to  the  experimental  data  than  the  Leitner’s  results  [9]. 

There  is  no  study  on  the  diffraction  patterns  due  to  a  point  source  incidence  although 
an  exact  solution  exists.  The  solution  is  expanded  in  terms  of  oblate  spheroidal  wave 
functions  which  are  not  simple  to  calculate.  Since  the  oblate  spheroidal  wave  functions  are 
not  available  in  any  mathematical  software  package,  a  computer  code  is  written  to  calculate 
the  functions  which  are  discussed  in  Chapter  4. 
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Pressure,  linear  scale 
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Total  pressure  field  amplitude  on  the  surface  of  the  disk 

ka  =  3 


Figure  3.20  The  BOR  simulation  of  the  disk  3 


Total  pressure  field  amplitude  on  the  surface  of  the  disk 

ka  =  4 


Figure  3.21  The  BOR  simulation  of  the  disk  4 
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Pressure,  linear  scale 


Total  pressure  amplitude  on  the  surface  of  the  disk 

ka  =  5 


Figure  3.22  The  BOR  simulation  of  the  disk  5 
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CHAPTER  4 

SPHEROIDAL  WAVE  FUNCTIONS 


Chapter  3  describes  BOR  in  detail  and  presents  the  BOR  simulation  results  compared 
with  the  experimental,  analytic  and  numerical  solutions  from  [5,  6,  7,  8,  9,  and  10]. 
However  those  results  are  valid  only  for  an  incident  plane  wave,  and  no  solution  is  available 
for  a  point  source  excitation.  In  order  to  verify  the  BOR  scattered  patterns  for  a  point  source 
excitation,  an  analytic  solution  for  an  axisymmetric  geometry  must  be  sought.  Studying  the 
oblate  spheroidal  geometry  is  a  good  way  to  verify  the  BOR  codes.  A  “fat”  spheroid 
represents  a  sphere,  and  a  thin  spheroid  becomes  a  disk. 

4.1  The  Oblate  Spheroidal  Geometry 

The  oblate  spheroidal  coordinates  (^,T|,<t))  shown  in  Figure  4.1  are  related  to  the 
rectangular  Cartesian  coordinates  (x,y,z)  by  the  transformation 

X  =  L d  +  l)(l  -  Tj^)  costj) 
y  =  ^d  l)(l  -Tl^)  sine!) 

z  =  ^d^Tl 

where  0<^<<»,  -1<ti<  1,  and  0  <  <()<  27C.  The  z-axis  is  the  axis  of  symmetry,  and  the 
interfocal  distance,  minor  axis,  and  major  axis  are  d,  d^  and  dV^^+l  ,  respectively  [6]. 
Because  of  the  axisymmetry  there  is  no  (j)  variation,  so  (})  is  assumed  to  be  zero  in  the  rest  of 
this  chapter. 

The  exact  solution  for  an  acoustically  hard  oblate  spheroid  due  to  a  point  source 
excitation  on  z-axis  (T|o  =  0)  is  computed  by  Bowman  et  al.  [7]  as 
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The  above  notations  are  consistent  with  the  Flammer’s  notations  [11].  The  coordinate  of  the 
oblate  spheroid  surface,  the  source  point,  and  the  field  point  are  denoted  as 
(^o>0,0),  and  respectively.  Also  and  take  the  minimum  and  maximum  values 

between  the  source  point  and  the  field  point  coordinates  respectively. 


Figure  4.1  The  oblate  spheroidal  coordinate  system 

The  oblate  spheroidal  coordinate  system  is  one  of  the  coordinate  systems  in  which  the 
scalar  wave  equation 

(vSk2)\|f  =  0  (4.3) 

is  separable  [1 1].  To  express  this  equation  in  spheroidal  coordinates,  the  following  relation 
for  metrical  coefficients,  h^,  hr),  and  h^  are  defined  by 
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dx2  +  dy2  +  dz2  =  d^^  +  dT)^  +  d(t)^ 


(4.4) 


where  the  scale  factors  are  defined  in  [12]  as, 


1/2 

h  -d. 

2\ 

4  +Ti2 

j 

.  h^-2 

i  l--n2  j 

[(l  l) 


1/2 


(4.5) 


With  the  use  of  the  expression  for  the  Laplacian  operator  in  orthogonal  curvilinear 
coordinates,  the  following  differential  equation  is  obtained  [12]. 

.2 


9ti 

=  I 

/ 

is  found  as 
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94  (4^  + i)(i -Ti2)a(t)^ 


■  +  C- 


\j/= 0 


(4.6) 


where  c  =  ^  kd.  By  the  usual  procedure  of  the  separation  of  variables,  the  solution  of  (4.6) 


¥mn  =  Smn(-ic,  Tj)  Rmn(-ic,i4)  exp(im(t)) 


(4.7) 


The  angle  and  radial  functions,  Sn,n(-ic,Ti)  and  Rnin(-ic,i4),  satisfy  the  ordinary  differential 
equations  [11] 


ati 


fe^+i)i-s„„{.ic,i4) 

fW-ic)-c2  4^  5"  1 

35  J 

h+K 

(1  -  ^  Smn(-iC,Tl)  +  ^n(-ic)  +  c2  ri^  - 

Ofi  J  j  .  .r|2 


^mn(”iC)i4)  ~  (^ 

Smn(“iC5”n)  =  0 


(4.8) 

(4.9) 


The  angle  and  radial  functions  are  the  key  functions  in  computing  (4.2).  All  of  the  variables 
in  (4.2)  are  evaluated  at  m  =  0,  and  the  rest  of  this  chapter  assumes  m  =  0.  Flammer 
describes  the  derivations  of  the  angle  and  radial  functions  in  detail  [11],  and  the  following 
two  sections  briefly  explain  his  derivations  since  understanding  the  derivations  are  crucial  in 
developing  a  computer  code. 
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4.2  The  Angle  Functions 

In  the  differential  equation  (4.9),  Son  is  obtained  from  the  expansions  in  the 
associated  Legendre  functions  of  the  first  kind  [11]: 

oo 

Son(-ic,ri)=  I'  d?"(-ic)  Pr(Tl)  (4.10) 

r=0,l 


The  prime  over  the  summation  sign  indicates  that  the  summation  is  over  only  even  values  of 
r  when  n  is  even,  and  over  only  odd  values  of  r  when  n  is  odd.  Substituting  (4.10)  in  (4.9) 
and  the  use  of  the  associated  Legendre  differential  equation  and  the  recursion  formulas  for 
the  associated  Legendre  functions  yield  the  following  recursion  formula  for  the  coefficients 
d?"  [11] : 


(r  +  2)(r+  l)c2  ^on 
(2r  +  3)  (2r  +  5) 


r(r+l)-Xon- 


r(r+l)-l  ^7 
(2r  -  1)  (2r  +3)  . 


d?"  + 


r  (r-1)  c^ 
(2r-2)  (2r-l) 


(r>0) 


(4.11) 


Finding  the  solution  to  the  above  recursion  formula  is  the  key  step  in  obtaining  the  angle 
functions  since  the  associated  Legendre  functions  are  relatively  well-known  functions.  A 
FORTRAN  subroutine  from  Numerical  Recipes  was  used  for  the  Legendre  functions  (see 

Appendix  C). 

For  convenience,  let 


Yr  =  r(r+l)-^c2 


1  -H 


1 

(2r  - 1)  (2r  +  3). 


(r^O) 


(4.12) 


r2  (r  - 1)^  c^ 

(2r  - 1)2  (2r  -  3)  (2r  +  1)  ’ 
_  r  (r  - 1)  c2  dp" 

■'(2r-l)(2r+l)d0n  ’ 


(r>2) 


(r>2) 


Substituting  above  notations  in  (4.1 1)  yields, 


(4.13) 

(4.14) 


Nr  = - - -  ,  (r>2)  (4.15) 

Yr  *  ^On  ■  Ni-+2 

and  reciprocally 

Nr+2  =  ^On  -  Yr  -  ^ ,  (r  >  2)  (4. 16) 

with  N2  =  lon  -  Yo,  N2  =  Xon  -  Yl- 

The  recursion  formula  (4.11)  must  be  convergent,  and  this  condition  leads  to  the 
iterative  computations  of  Nj  from  (4.14)  for  large  values  of  r.  Conversely  from  the  initial 
values  N2  and  N3,  Nf  is  calculated  from  (4.15)  for  small  values  of  r.  When  c  is  small,  the 
dominant  coefficient  for  given  n  is  d^"  [11].  Thus  it  is  convenient  to  obtain  the  ratios  of  the 
coefficients  d?"/dS".  It  is  observed  that  the  series  of  the  coefficient  ratios  obtained  from 
(4.15)  is  accurate  only  for  large  r,  and  the  series  from  (4.14)  is  accurate  only  for  small  r. 
Therefore  for  r  less  than  n,  the  coefficient  ratios  from  (4. 14)  are  taken,  and  for  r  greater  than 
n,  the  coefficient  ratios  from  (4.15)  are  taken. 

The  calculation  of  the  coefficients  are  based  on  the  accurate  value  of  ^on-  The 
eigenvalue  for  a  small  c^  is  obtained  from  a  series  in  powers  of  c^  [11]. 

00 

?lOn(-ic)=  I  (-l)MSc21c  (4.17) 

k=0 

where  ig"  =  n  (n+1) 

lOn  =  i[l  + - 1 - 

^  2  [  (2n  -1)  (2n  +3). 

jOn  _ _ -  (n  +1)^  (n  +  2)2 _ ^ _ n2  (n  -  1)2 _ 

2  (2n  -  1)  (2n  +  2)  (2n  +  3)3  (2n  +  5)  2  (2n  -  3)  (2n  -  1)3  (2n  +  1) 

|0n^ _ n2  (n  - 1)^ _ ^ _ (n  +1)^  (n  +  2)^ _ 

^  (2n  -  5)  (2n  -  3)  (2n  - 1)3  (2n  +  1)  (2n  +  3)  (2n  - 1)  (2n  +  2)  (2n  +  3)3  (2n  +  5)  (2n  +  7) 

(n  - 1)2 _ (n  +  1)2  (n  +  2)2 _ 

(2n  -  5)2  (2n  -  3)  (2n  -  1)^  (2n  +  1)  (2n  +  3)^  (2n  - 1)2  (2n  +  1)  (2n  +  3)^  (2n  +  5)  (2n  +  7)^ 
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(2n-  1)^  (2n-2)^  (2n  -  3)^ _ (n  +  1)^  (n  +  if  (n  +  3)^  (n  +  4)^ 

(2n  -  7)  (2n  -  5)^  (2n  -  3)3  (2n  -  1)^  (2n  +  1)  (2n  +  1)  (2n  +  3)^  (2n  +  5)^  (2n  +  if  (2n  +  9) 

(n  +  1)^  (n  +  If _ n4  (n  - 1)^ 

”  (2n  +  1)2  (2n  +  3)^  (2n  +  5)2  (2n-  -  3)2  (2n  -  1)^  (2n  +  1)2 


_ (n  -  1)2  (n  +  1)2  (n  +  2)2 _ 

(2n  -  3)  (2n  - 1)^  (2n  +  1)2  (2n  +  3)^  (2n  +  5) 


(418) 


In  the  numerical  computation,  truncating  the  summation  up  to  k  =  4  yields  sufficiently 
accurate  eigenvalues  for  small  c^.  For  a  large  argument,  an  asymptotic  expansion  is  used 
[11]. 


?lOn(-ic)  =  -  c2  +  2c  (2v  +  1)  -  2v  (v  +  1)  -1  +  Aon 

where  v  =  —  for  n  even  and  v  =  for  n  odd.  Aon  is  defined  as 
2  2 


(4-19) 


A«„  =  S  P?  c-k 

k=l 


P?"  =  -  2-3  q  (q2  +  1)  ,  P?"  =  -  2-6  [5q4  +  10q2  +  l] 

P?  =  -  2-9  q  [33q4  +  1 14q2  +  3?]  ,  p?  =  -  [63q6  +  340q4  +  239q2  +  14]  (4.20) 


where  q  =  n  +  1  for  n  even  and  q  =  n  for  n  odd. 

The  power  series  and  asymptotic  expansion  methods  are  not  sufficient  for 
intermediate  values  of  c2,  and  the  eigenvalue  should  be  refined  by  Boukamp’s  method  of 
approximation.  The  eigenvalue  is  obtained  from  a  transcendental  equation  U  which  is 
originated  from  equations  (4.15)  and  (4.16)  [1 1] 

U(^n)  =  Ul(^n)  U2(^n)  =  0  (4-21) 

where 
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U,(?iOn)  =  -%  +  hn - - 

Yn-2  "  ^On  ■  N[|_2 


and  U 1  (Xon)  = - - 

Yn+2  ■  ^On  ■  Nn+4 


Note  that  Nn.2  is  evaluated  by  the  equation  (4.16)  and  Nn+4  is  evaluated  by  (4.15).  Let  X.on 

be  the  approximate  eigenvalue  calculated  from  either  of  the  above  two  methods,  and  5^on 

the  difference  between  the  actual  eigenvalue  and  the  approximate  value.  | 

0  =  U(kon)  =  U(X,i)‘„^  +  6?l0n)  »  U(XE)n^  +  5U  (4.22) 

i 

! 

by  finding  the  first  variation  of  U  due  to  the  variation  6Xon-  If  the  variation  5A.on  is  made  in  i 

the  eigenvalue,  the  variation  on  (4.16)  is 


5N„2  =  8Xo„  +  ;%8N,  (4.23) 

(N,)2 

Likewise,  the  variation  on  (4.15)  is 

5Nr  =  ^  [5Nr+2  -  Skori  (4.24) 

Pr 


By  iteration  of  (4.23)  and  (4.24)  the  variations  of  U  are  obtained 

Pn  Pn-2  Pn-4 


5Ui  =  6?l0n 


5U2  =  8>.0n 


1  +  P"  ^  Pn  Pn-2  ^ 


{N„)2  (Nn)2(N„.2)2  (Nn)2(N„.2)2(N„.4)2 

(Nn)^  ,  (Nn)^(Nn^2)^  ,  (Ng)^  (Nn^2f  (Nn^4F  , 
Pn  Pn  Pn+2  Pn  Pn+2  Pn+4 


(4.25) 


(4.26) 


Substituting  (4.25)  and  (4.26)  in  (4.22),  5A,on  is  found  to  be 

-Ui(?.y„^)-U2(>.d 


8X,0n  ' 


1  I  Pn  I  PnPn-2  ^ 

(Nn)2  (Nn)2(Nn-2)2 


(NniZ  +  lNoiafitOi 


L  P 


n+2 


Pn+2  Pn+4 


(4.27) 


The  new  eigenvalues  obtained  by  the  above  method  yield  exactly  same  eigenvalues  listed  in 
[12].  Thus  the  Boukamp’s  method  leads  to  remarkably  accurate  eigenvalues,  and  this 


i 


I 
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accuracy  is  the  key  in  calculating  the  coefficients. 

Once  the  coefficient  ratios  are  determined  accurately,  the  actual  values  of  the 
coefficients  can  be  obtained  in  terms  of  an  arbitrary  coefficient  value.  Flammer  uses  a 
normalization  scheme  [12]  such  that  each  spheroidal  angle  function  reduces  exactly  to  the 
corresponding  associated  Legendre  function  when  c  becomes  zero.  This  normalization 
follows  Chu  and  Stratton’s  normalization  scheme  [12]  except  that  the  normalization  is  earned 
out  at  T|  =  0.  Thus  using  the  normalization  relations  below,  the  coefficients  are  completely 

determined. 

y'HiilLdpn  for  n  even  (4.28) 


yi  (-1)  2  (r  +  i)!  = — ('^)  ^  for  n  odd 


(4.29) 


Once  the  expansion  coefficients  are  determined,  Nnn  from  (4.2)  is  easily  found  as  below 

[11]. 


oo 

Non  =  2  I' 


r=0,l 


r!  {d?"f 
(2r+  l)r! 


(4.30) 


4.3  The  Radial  Functions 

The  radial  functions  satisfy  the  differential  equation  (4.9).  The  eigenvalue  in  (4.9)  is 
identical  as  in  (4.8)  of  the  angle  function.  The  radial  functions  are  found  in  terms  of  the 
spherical  Bessel,  Neumann,  and  Hankel  functions  [8]. 


I' 


Lr=0.1 


1-1 


X*  P"d?"Z^P\c4) 

r=0,l 


(4.31) 


where 


Z<P\z)  =  y^Jr+l(z)  forp=l 

=  forp  =  2 


(4.32) 
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J  and  Y  are  Bessel  and  Neumann  functions,  respectively.  Similar  to  Hankel  functions,  the 
radial  functions  are  related  likewise  [12] 

R®(-ic,iy  =  C(-ic,i«  +  i  RS(-ic,i|) 

(4.33) 

RS;;(-ic,i^) = Ri‘„Vic,i^)  -  i  Rg;(-ic,i^) 

(4.33)  makes  defining  the  radial  functions  easier.  The  derivative  of  Ro[,(-ic,  i^)  is  easily 
obtained  by  taking  the  first  derivatives  of  RQj|(-ic,i^)  and  RQ^jj\-ic,i^) . 

R[jJ'(-ic,i4)  = 

However  calculating  the  radial  function  of  the  second  kind  for  small  values  of  %  is  not 
suitable  since  the  Neumann  function  converges  poorly  for  small  arguments.  For  small  values 
of  a  better  numerical  method  should  be  used  to  calculate  the  radial  function  of  the  second 
kind. 

The  angle  functions  and  the  radial  functions  are  related  by  joining  factors  Kq^j'  [12], 
S[f„\-ic,i^)  =  Kg;\-ic)  RS(-ic,i^)  for  p  =  1,  2  (4.35) 

where 


I'  d?"  I'  d?"  ZSp^'(c4) 

r=0,l  J  r=0,l 


(4.34) 
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Based  on  the  above  relations,  Flammer  has  derived  power  series  expansions  for  the  oblate 
radial  functions  as  below  [11]. 


R®(-ic,  i^)  =  Qon(-ic)  Ro„Vic,i^)  I  tan'k^)  -  f  +  gOn(-ic,i^) 


2J 


(4.37) 


with 


Qon(-fc)  =  ag°(-ic)  for  n  even 


0- 

QonOO-  c  '-'''-',>[2-r(-r)!p 


On/  _Ji — 2r)^ —  ^ 


a^C-ic)- 


where  ao"(-ic)  =  ^  ^ •  The  coefficients  c^U  are 

[c°of 


(k!)2  r=k 


4k  =  7^1  (-4  (r  + 1)  d2r+i  for  n  odd 

(k!)2  r=k  ^  2/k 


(4.38) 


(4.39) 


where  (r)k  =  r  (r  +  1)  (r  +  2)  •  •  •  (r  +  k  - 1),  (r)o  =  1 . 

The  function  gon  in  (4.37)  satisfies  the  inhomogeneous  radial  equation  [11] 

-d-  +  l)  -  Xon(-ic)  +  c^^^  gOn(-ic,i^)  =  -  2  Qon(-ic)  RQ^\-ic,i^) 

d^  d^  J 


(4.40) 


It  is  convenient  to  expand  gon  in  the  form 


g0n(-ic,iO  =  I  b5?  for  n  even 

r=0 

gOn(-ic,i^)  =  I  'C'  for  n  odd 

r=n 


(4.41) 


since  the  power  series  expansions  are  used  for  the  values  of  ^  which  are  close  to  zero  [1 1]. 
After  substituting  (4.41)  in  (4.40)  and  rearranging,  the  following  recursion  formulas  for 
are  obtained  [1 1]. 
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(4.42) 


(2r  +  2)  (2r  +  3)  +  [(2r  +  1)  (2r  +  2)  -  Xo„(-ic)]  BJ  +  c2  b5,".2 

oo 

=  -  2  Qon(-ic)  [Kon^-ic)]  '  S  c^k(-ic)  2k  +  1)  for  n  even, 

k=r+l 


(2r  +  1)  (2r  +  2)  b5?,2  +  [2r  (2r  +  1)  -  X<,„(-ic)]  Bg  +  c2  B§” , 


(4.43) 


=  -2Qo„(-ic)[i-'K^'„>(.ic)]'' 


Ic^E(-ic)(2k+l)©-  2  c§E(-ic)2k(*‘;>) 

k=r  k=r+l 


for  n  odd. 


The  symbol  (^)  denotes  the  binomial  coefficient  ■■  . .  The  initial  coefficients  Bn" 

‘  (k-r)!r! 


are 


found  in  [11]  as 

Bo"  =  [c  Ron(-ic40)]  ^  -  Qjn(-ic)  Rjj„\-ic,i0)  for  n  even 
Bq"  =  [c  Rq^^  (-ic,iO)]  ^  for  n  odd 

where 


(4.44) 


Ro„\-ic,iO)  =  --.^0  ^  for  n  even,  RQ„\-ic,iO)  =  0  for  n  odd 

oo 

I'd?"(-ic) 

r=0 


Rj)^„)’(-ic,iO)  =  ^ 


n-1 


dr(-ic) 


3X'  d?"(-ic) 

r=0 


for  n  odd. 


^On  =  0  for  n  even 


(4.45) 


The  remaining  coefficients  are  completely  determined  from  the  recursion  formulas  in  (4.42) 
and  (4.43). 

The  derivative  of  Rq^^^^  is  determined  by  simply  taking  a  derivative  of  (4.37). 


=  Qon(-ic)  Ron  [tan-^(^)  - 


4‘^"(-‘c)5^'‘°-‘^)4g'0n(-ic.i4)(4.46) 


E2  1 
^  + 1 


where  g'on  is 
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oo 


I(2r  + 

r=0 


for  n  even 


g'«„(-ic,©=l2rB5?5''-'  for  n  odd 


r=n 


(4.47) 


Therefore  the  numerical  means  of  calculating  the  radial  functions  are  completely  determined. 


4.4  Analytic  Solutions 

The  exact  solution  in  (4.2)  have  been  calculated  in  terms  of  the  angle  and  radial  oblate 
spheroidal  functions.  Although  (4.2)  is  a  short  equation,  writing  a  FORTRAN  code  for  it  has 
been  challenging.  There  is  a  limited  listing  of  the  eigenvalues,  coefficients,  angle  functions 
and  radial  functions  in  [11]  and  [12].  A  portion  of  those  values  listed  in  [1 1]  and  [12]  have 
been  verified  when  debugging  the  code.  In  spite  of  the  inadequate  information,  the  code 
calculates  satisfactorily.  Figures  4.2  and  4.3  compare  the  analytic  solutions  obtained  by  the 
code  and  the  corresponding  BOR  simulation  results  for  thin  disks.  The  exact  solution  in 
(4.2)  involves  a  summation  from  zero  to  infinity  for  n,  but  a  summation  from  zero  to  30 
proves  to  be  adequate  for  the  code.  Comparing  the  analytic  solution  has  proved  that  the  BOR 
code  is  valid  for  a  point  source  excitation  as  well. 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  4.2  The  BOR  simulation  and  analytic  diffraction  solution  of  the  disk  1 


Total  pressure  field  amplitude  in  the  azimuthal  plane 


Figure  4.3  The  BOR  simulation  and  analytic  diffraction  solution  of  the  disk  2 
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CHAPTER  5 
CONCLUSIONS 


In  this  thesis  numerical  methods  for  calculating  diffractions  of  acoustic  waves  have 
been  investigated.  Diffraction  patterns  have  been  accurately  calculated  using  the  MOM  code 
for  simple  objects  small  compared  to  their  wavelength.  Chapter  2  introduces  the  MOM  and 
presents  the  simulation  results.  The  MOM  code  is  most  useful  for  studying  diffractions  of 
rectilinear  objects,  and  it  is  not  suitable  for  analyzing  curved  objects  due  to  the  unfriendly 
meshing  scheme.  The  BOR  code  is  written  to  compensate  for  the  shortcomings  of  the  MOM 
code.  Diffraction  patterns  of  axisymmetric  objects  with  axisymmetric  boundary  conditions 
can  be  efficiently  calculated  by  the  BOR  code.  The  diffraction  patterns  of  spheres,  parabolic 
reflectors,  and  thin  disks  have  been  successfully  calculated,  and  the  simulation  results  are 
shown  in  Chapter  3.  Diffractions  due  to  a  plane  wave  incidence  are  easily  compared  with 
existing  references.  For  a  point  source  incidence,  no  graphical  data  are  available,  and  only 
the  analytic  solution  is  expressed  in  terms  of  the  spheroidal  wave  functions.  In  Chapter  4  the 
oblate  spheroidal  functions  are  discussed  in  detail.  A  code  is  written  to  calculate  the  analytic 
solution,  and  the  calculation  results  have  verified  the  validity  of  the  BOR  code  for  a  point 
source  incidence. 

Note  that  the  boundary  conditions  must  also  be  axisymmetric  for  the  BOR 
simulations.  All  of  the  BOR  simulation  results  in  this  thesis  are  done  in  the  full  space.  For  a 
plane  wave  incidence,  the  direction  of  the  propagation  is  parallel  to  the  axis  of  symmetry,  and 
for  a  point  source  incidence,  the  source  is  located  on  the  axis.  Analyzing  objects  in  a  half 
space,  such  as  spheres  above  the  flat  ground,  can  also  be  done  using  the  BOR  code.  Placing 
acoustic  images  of  the  original  objects  in  the  full  space  is  equivalent  to  the  original  objects  in 
the  half  space  which  loses  the  symmetry. 

The  BOR  code  can  further  be  improved  if  it  is  modified  to  account  for  off-axis 
incidence.  The  current  code  assumes  that  there  is  no  field  variation  on  a  plane  perpendicular 
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to  the  symmetric  axis.  For  example,  when  the  plane  wave  is  incident  at  an  angle  from  the 
axis  of  the  symmetry,  the  incident  fields  on  the  surface  of  the  object  will  vary  accordingly. 
The  scattered  fields  will  have  the  same  variation  due  to  the  variation  of  the  incident  fields. 
Finding  out  of  a  way  to  incorporate  this  variation  in  the  BOR  code  will  broaden  its  usage  in 
analyzing  the  diffractions. 


54 


REFERENCES 


[1] 

[2] 


C.A.  Balanis,  Advanced  Engineering  Electromagnetics,  Wiley,  New  York,  1989. 


G  W  Swenson  Jr.,  E.R.  Sandeen,  L.L.  Pater,  and  H.C.  Zhuang,  "The  potential 
for  mitigation’ of  gun  blast  noise  through  sheltering  of  the  source,"  Technical  Report 
N-92/09,  United  States  Army  Corps  of  Engineers,  Construction  Engmeenng 
Research  Laboratory,  April,  1992. 


[3]  N.  Morita,  N.  Kumagai,  and  J.  R.  Mautz,  Integral  Equation  Methods  For 
Electromagnetics,  Artech  House,  Boston,  1990. 


[4] 

[5] 

[6] 

[7] 


J.  M.  Jin,  The  Finite  Element  Method  In  Electromagnetics,  Wiley,  New  York, 
’ 1993.  ’ 

J  W  Benson  Y.  L.  Li,  and  G.  W.  Swenson,  Jr.,  “A  baffle-type  directional 
microphone,”  J.  Acoust.  Soc.  Am.,  vol.  95,  pp.  2536-2538,  1994. 


A.  F.  Seybert,  B.  Soenarko,  F.J.  Rizzo,  and  D.J.  Shippy,  “A  special  integral  equation 
formulation  for  acoustic  radiation  and  scattering 

boundary  conditions,”  J.  Acoust.  Soc.  Am.  ,  vol.  80,  pp.  1241-1247,  1986. 

J.  J.  Bowman,  T.B.A.  Senior,  P.L.E.  Uslenghi,  Electromagnetic  and  Acoustic 
Scattering  by  Simple  Shapes,  Hemisphere  Publishing  Corp.,  New  York,  1987. 


rSlS  Wahlstrom,  “The  Parabolic  Reflector  as  an  Acoustic  Amplifier,”  J.  Audio  Eng.  Soc., 
vol.  33,  pp.  418-429,  1985. 

[9]  A  Leitner,  “Diffraction  of  Sound  by  a  Circular  Disk,”  J.  Acoust.  Soc.  Am. ,  vol.  21, 
pp.  331-334,  1949. 

riO]  F.  M.  Wiener,  “The  Diffraction  of  Sound  by  Rigid  Disks  and  Rigid  Square  Plates,”  J. 
Acoust.  Soc.  Am.,  vol.  21,  pp.  334-347,  1949. 


[11]  C.  Flammer,  Spheroidal  Wave  Functions,  Stanford  University  Press,  Stanford, 
1957. 


[12] 


M.  Abramowitz  and  I.  A.  Stegun,  Handbook  Of  Mathematical  Functions  With 
Formulas,  Graphs,  And  Mathematical  Tables,  Dover  Publications,  Inc.,  New 
York,  1970. 


55 


USACERL  DISTRIBUTION 


Chief  of  Engineers 

ATTN:  CEHEC-IM-LH  (2) 

ATTN:  CEHEC-IM-LP  (2) 

ATTN:  CEMP 
ATTN:  CEMP-CE 
ATTN:  CEMP-EA  (2) 

ATTN:  CEMP-ZM 
ATTN:  CERD-L 
ATTN:  CERD-M  (2) 

ATTN:  CERD-ZA 

ACS(IM)  22060 

ATTN:  DAIM-FDP 

CECPW  22310-3862 
ATTN:  CECPW-E 
ATTN:  CECPW-FT 
ATTN:  CECPW-ZC 

Norwegian  Defense 

ATTN:  Chief  Test  &  Dvlpmt  Section 

HQ  USAREUR  &  7th  Army 
ATTN:  AEAEN-EH 
ATTN:  Unit  29351 

US  Military  Academy 
ATTN:  MAEN-A 
ATTN:  Civil  Div  Director 
ATTN:  Dept  of  Geo  &  Env  Engr 

US  Army  Combat  Sys  Test  Activity 
ATTN:  STECS-SO-R  21010-5059 

Commander  FORSCOM 

ATTN:  FCEN-RDF  30330-6000 

CEWES  39180 
ATTN:  Library 

CECRL  03755 
ATTN:  Library 


US  Army  ARDEC  07806 
ATTN:  SMCAR-ISE 

Engr  Societies  Library  10017 
ATTN:  Acquisitions 

US  Army  Environmental  Center 
ATTN:  SFIM-AEC-NR  21010 
ATTN:  SFIM-AEC-CR64152 
ATTN:  SFIM-AEC-SR  30335-6801 
ATTN:  AFIM-AEC-WR  80022-2108 

National  Guard  Bureau  203 10 
ATTN:  NGB-ARI 

US  Military  Academy  10996 
ATTN:  MAEN-A 
ATTN:  Facilities  Engineer 

Naval  Facilities  Engr  Command 
ATTN:  Facilities  Engr  Command 
Code  20YAZ  (2) 

US  Army  CHPPM 

ATTN:  MCHB-DC-EEN  (5) 

US  Gov’t  Printing  Office  20401 
ATTN:  Rec  Sec;Deposit  Sec  (2) 

Nat’l  Institute  of  Standards  &  Tech 
ATTN;  Library  20899 

Defense  Construction  Supply  Center 
ATTN:  DCSC-WI 43216-500 

Defense  Tech  Info  Center  22060-6218 
ATTN:  DTIC-0  (2) 

49  +  12 
5/97 


56 


